Lattice Theory of Pseudospin Ferromagnetism in Bilayer Graphene: 
Competing Orders and Interaction Induced Quantum Hall States 

Jeil JungQ Fan Zhang, and Allan H. MacDonald 

Department of Physics, University of Texas at Austin, USA 
(Dated: November 30, 2010) 

In mean-field-theory bilayer graphene's massive Dirac fermion model has a family of broken inversion sym- 
metry ground states with charge gaps and flavor dependent spontaneous inter layer charge transfers. We use a 
lattice Hartree-Fock model to explore some of the physics which controls whether or not this type of broken 
symmetry state, which can be viewed as a pseudospin ferromagnet, occurs in nature. We find that inversion 
symmetry is still broken in the lattice model and estimate that transferred areal densities are ~ 10~ 5 electrons 
per carbon atom, that the associated energy gaps are ~ 10~ 2 eV, that the ordering condensation energies are 
per carbon atom, and that the energy differences between competing orders at the neutrality point are 
~ 10 _9 eV per carbon atom. We explore the quantum phase transitions induced by external magnetic fields and 
by externally controlled electric potential differences between the layers. We find, in particular, that in an ex- 
ternal magnetic field coupling to spontaneous orbital moments favors broken time-reversal-symmetry states that 
have spontaneous quantized anomalous Hall effects. Our theory predicts a non monotonic behavior of the band 
gap at neutrality as a function of interlayer potential difference in qualitative agreement with recent experiments. 
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I. INTRODUCTION 

Recent experimental progress^ in isolating and measur- 
ing the electronic properties of graphene single and multilay- 
ers has opened a new topic in two-dimensional electron sys- 
tem (2DES) physics 3 . Electronic wavefunctions in graphene 
systems are often described using a pseudospin language in 
which the spinors specify wavefunction components on differ- 
ent sublattices. Although the properties of graphene 2DES's 
can often be successfully described using an effective non- 
interacting electron model, studies of electron-electron in- 
teractions effects have revealed some qualitative differences 
compared to ordinary 2DES's 4,5 that are related to these sub- 
lattice pseudospin degrees-of-freedom. 

In the case of AB stacked bilayer graphene, there are four C 
sites and four 7C-orbitals per unit cell, but two of these are re- 
pelled from the neutral system Fermi level by interlayer hop- 
ping. This circumstance leads to a low-energy massive chiral 
fermion model 6 with two-component spinors, and a crystal- 
momentum p dependent pseudo-magnetic field with a magni- 
tude that varies as p 2 and an orientation angle twice as large 
as the momentum orientation angle §p. Neutral system states 
have one occupied pseudospin for each distinct set of momen- 
tum, spin, and valley labels. Recently 7 Min et al. pointed 
out that when Coulombic electron-electron interactions are 
added to the massive chiral fermion model, the mean-field the- 
ory ground state pseudospins of each spin-valley flavor break 
symmetry by rotating out of the x — y plane, developing z 
components with a common spontaneously chosen sign and 
magnitudes which are larger at small p. Because the two- 
sublattices from which the two-band-model pseudospins are 
constructed are located in opposite layers, the broken symme- 
try transfers charge between layers. It is therefore character- 
ized in momentum space by a vortex with vorticity v = 2 and 
a flavor-dependent core polarized along one of the polar di- 
rections, and in real space by flavor-dependent uniform layer 
polarization. Different members of the family of states are dis- 



tinguished by the spin-valley flavor dependence of the sense 
of layer pseudospin orientation in the momentum-space vor- 
tex cores. The origin of the broken symmetry, which we refer 
to here as pseudospin ferromagnetism, 7 is the p 2 pseudospin- 
splitting at small p, which leads to infrared divergences^ in 
particle-hole polarization loops, combined with the frustrat- 
ing effect of pseudospin chirality which leads to relatively 
stronger exchange interactions for £-polarized pseudospins. 

The broken symmetry states can be classified*^ in a man- 
ner which highlights their valley and spin dependent mo- 
mentum space Berry curvatures, which are in turn closely 
related 1112 to the Hall responses of the system. Indeed, re- 
cent experiments 13 demonstrate that bilayer graphene exhibits 
a quantized quantum Hall effect in the absence of an exter- 
nal magnetic field. Over and above the basic science inter- 
est deriving from this presently unique property, it is possible 
that pseudospin order in bilayer graphene could play a role in 
graphene electronics 14 by introducing hysteresis and dramati- 
cally enhancing the influence of gates on conductance. 

In this article we report on a study of pseudospin ferromag- 
netism in a jr-orbital tight-binding model for bilayer graphene, 
which is conveniently able to capture and establish the role 
of some bilayer graphene 7C-band features neglected in the 
massive chiral fermion model. We estimate that the density 
shift for each flavor is ~ 10~ 5 electrons per carbon atom, 
that the gaps are ~ 10~ 2 eV, that the total condensation en- 
ergy is ~ 10~ 7 eV per carbon atom, and that the energy dif- 
ferences between competing ordered states is ~ 10~ 9 eV per 
carbon atom. The energy differences between competing or- 
dered states, which are between one to two orders or mag- 
nitude smaller than the condensation energy, are sensitive to 
roughly estimated lattice scale details of the model we em- 
ploy. The competing states are classified in the first place as 
either anomalous Hall states, in which opposite valleys are po- 
larized toward opposite layers, or valley Hall states in which 
they are polarized toward the same layer. We find that val- 
ley Hall states have slightly lower energies. Because anoma- 
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lous Hall states have spontaneous orbital magnetism they are 
favored by external magnetic fields. We estimate that, be- 
cause the energy differences between states are very small, 
extremely weak fields are sufficient to induced phase transi- 
tions between valley Hall and anomalous Hall states. On the 
other hand potential differences between layers favor valley 
Hall states. When the spin degree of freedom is accounted for 
we find that the spontaneous energy gap in neutral bilayers 
first decreases and then increases with potential difference, in 
qualitative agreement with experiment. 

In the following two sections we explain the model Hamil- 
tonian we employ and some of the technical details of the 
calculations we have carried out. The main results are pre- 
sented in the Section IV where we introduce the topological 
classification of solutions and use it to discuss the properties 
of broken symmetry solutions for vanishing, moderate, and 
strong externally controlled electric potential differences be- 
tween the layers. (Below we refer to interlayer potential dif- 
ferences, which always play a key role in graphene bilayer 
physics, as potential biases.) We also discuss coupling to an 
external magnetic field which can drive transitions between 
competing states, some of which have broken time-reversal- 
symmetry with associated anomalous Hall effects and orbital 
magnetism. 

We close the paper with a summary and some suggestions 
for further research. 



II. FOUR-BAND ^-ORBITAL TIGHT-BINDING MODEL 
OF BILAYER GRAPHENE 

We describe bilayer graphene using a lattice model with 
one atomic 2p z orbital per carbon site. We write the model's 
Bloch basis states in the form 



Ykr(r) = 4eJ>*< B '- 



x *' (r — R; — Tjc) , 



(1) 



where N is the total number of unit cells in the system, <j> (r) 
is the band's Wannier wavefunction, and K labels the carbon 
site with position T K relative to a the triangular lattice vec- 
tor R,-. (We comment later on the possible role of screen- 
ing effects from the p and s orbitals which form the d and 
a* bonds neglected in this model.) Following the convention 
used in Refs. l6lfl5ll . we use the notations A, B, A, B for the 
four sublattice indexes K, where B and A are the opposite- 
layer near-neighbor-pair sites. With this convention, the four 
band tight-binding model Hamiltonian of a graphene bilayer 
is: 
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with a = 2A6A arises from a sum over the three near-neighbor 
hops within a layer. We have neglected differences in on- 
site energies and next nearest neighbor hopping processes 
which give rise to electron-hole asymmetry and do not play 
an important role in pseudospin ferromagnetism. The tight- 
binding model parameters Y should not be confused with 
the Slonczewski- Weiss, McClure 16 model parameters for bulk 
graphite, despite the obvious similarities in notation. In our 
calculations we adopt conventions similar to those of Ref. lfl7ll 
for bilayer graphene, taking the values 70 = — 3.12 eV, Yi = 
-0.377, 73 = -0.29 eV and 74 = -0.12 eV for the hopping 
parameters. Only the intralayer nearest neighbor (70) process 
and interlayer tunneling (71) process are retained in the min- 
imal tight-binding model. The trigonal warping (73) process 
which connects the A and B sites is responsible for the lead- 
ing circular symmetry breaking near the valley points, while 
the (74) process which connects A and A sites influences the 
intralayer charge imbalance between sublattices A and B. 



III. LATTICE MODEL MEAN-FIELD THEORY 

Because of the importance 5 of non-local exchange 
in graphene systems, a Hartree-Fock mean-field theory 
approximation 18 is a natural first step in considering electron- 
electron interaction effects. When Coulomb interactions are 
added to the 7T-band tight-binding model the interaction terms 
in the mean-field Hamiltonian take the form: 

V HF = £ uP'Nv&ch ~ W&^cw (4) 
kJU' 

where A is a composite label for sublattice K and spin a. The 
first term on the right hand side of Eq. (01 is the Hartree term: 
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where G is a reciprocal lattice vector. The second is the Fock 
(exchange) term: 
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In Eq.© and Eq.® the two-dimensional Coulomb interac- 
tion V KK (q) = 2ne L I (|q| e r ) when K and k' refer to the same 
layer and (2ne 2 / (|q|e r )) exp [ — jqj c] when k and k' refer to 
the opposite layers. Here e r is the relative dielectric constant, 
c = 3.35A is the interlayer separation, A is the total area of the 
graphene sheet, and we use 
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as a form factor which accounts for the spread of the ^-orbital 
charge on each site. This simple form assumes an isotropic 
site-localized charge distribution. Eq. © was obtained by 
Fourier transforming the radial charge distribution of a hydro- 
genic 2p orbital. The use of ro = oo = flo/V^O would yield 
a root mean square radius corresponding to the covalent ra- 
dius of the carbon atom oq = OJlA. If we consider screening 
from the a band electrons neglected in our model and the fact 
that the charge density distribution of a p z orbital is far from 
spherical we expect that larger values of ro, which effectively 
reduce onsite repulsion, would be more appropriate. For most 
of our calculations we have therefore used the value ro = 32o- 

As explained in the introduction, pseudospin ferromag- 
netism in bilayer graphene can be neatly described using the 
two-band massive chiral fermion model. This approach has 
two shortcomings which the present calculation is intended 
to alleviate. First of all, the model has to rely on a crude 
ultraviolet-cutoff to account for the limited range of energy 
~ Y\ over which it is applicable. At moderate interaction 
strengths the amount of charge transferred between layers de- 
termined in the massive chiral fermion model calculation is 
strongly influenced by this cut-off. Secondly, the calculation 
described in Ref. (0) relies on the model's circular symme- 
try for a number of simplifications. When direct hopping be- 
tween the low-energy A and B sites, the 73 process, is included 
in the Hamiltonian the model's Fermi lines are no-longer cir- 
cular and the continuum model loses some of its attractive 
simplicity . 19 ! 20 These processes are known to be essential at 
very weak interaction strengths since they remove the infrared 
divergences 8,21 responsible for the instabilities of the massive 
chiral fermion model. The present lattice model reduces to 
the continuum model at low energies, accounts naturally for 
the limited validity range of two-band models by retaining all 
four it -bands, and can deal with the loss of circular symmetry 
without any additional complication. 

The main challenges which arises in practical implementa- 
tion of the lattice model mean-field-theory lie in the numerical 
optimization of a problem in which the small portion of the 
Brillouin-zone close to one of the valley points plays a domi- 
nating role and must be sampled densely. It is essential that we 
have sufficiently dense fc-point sampling near the Dirac points, 
and at the same time sample the full Brillouin zone. The in- 
crease of computational load with fc-point sampling density 
is particularly rapid in the Hartree-Fock calculations because 
the Hartree-Fock matrix element at one fc-point depends on 
the occupied wavefunctions at all other fc-points. In Fig. Q] 
we illustrate the honeycomb lattice reciprocal space primitive 
cell and the fc-point sampling scheme we have chosen. In an 
effort to achieve a satisfactory compromise between compu- 
tational load and accuracy we, first of all, make use of the 
hexagonal symmetry inherent in the problem. This allows us 
to limit our calculations to 1 /6th of the total Brillouin zone 
area when we distinguish K and K' valleys, and 1/1 2th of the 
total area when we do not. In addition, instead of using a uni- 
form grid in the whole Brillouin zone we use a Appoint mesh 
which is denser near the Dirac point. In our calculations we 
have used 18 x 18, 32 x 32 or 42 x 42 coarse grids for the 
full primitive cell and multiplication factors of up to 128 for 
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FIG. 1: Left panel: Band structure of graphene represented in the 
primitive zone defined by the reciprocal lattice vectors b\ and i>2 of 
the honeycomb lattice. The equilateral triangle represents the region 
in the primitive cell we need to sample in order to fully describe a 
system in which the symmetry between K and K' valleys is broken. 
The smaller triangles enclose inequivalent regions in £-space asso- 
ciated with K and K' valleys. Right panel: An example of coarse 
sampling of fc-points in the primitive cell supplemented by a finer 
grid in hexagons generated around those coarse ^-points near the K 
and K 1 valleys. The illustrated example consists of 32 x 32 coarse 
points and a 16-fold enhancement of density resulting in an effective 
sampling of 512 x 512 points in the primitive cell near the valleys. 



the finer Appoint mesh near the valley centers. The fine mesh 
densities were either 2304 x 2304 or 2048 x 2048 for typi- 
cal calculations at zero potential bias, and 672 x 672 when 
smaller densities were enough to converge the calculations in 
the case of strongly biased bilayers. Dense Appoint meshes 
were normally employed within ~ 0.5/a of a Dirac point, 
and wider regions were used for strong bias cases. (Here 
a = 2A6A is the lattice constant of the triangular periodic lat- 
tice structure of graphene.) Nevertheless £-point sampling ap- 
proximations remain the main source of numerical inaccura- 
cies. The self-consistent-field calculations were iterated until 
convergence to 14 significant figures in the total energy per 
electron was achieved for a given choice of fc-point sampling. 
In order to resolve the small energy differences between so- 
lutions of the self-consistent field equations corresponding to 
different states we needed to compare results obtained with 
the same £-point sampling scheme. The £-point sampling we 
used achieved convergence to within a few parts per thousand 
for charge density differences and band gaps. 



IV. LATTICE MODEL PSEUDOSPIN FERROMAGNET 

The inversion symmetry breaking instability in bilayer 
graphene occurs nearly independently at the two points 
and, within mean-field-theory, entirely independently for 
each electron spin. 7 It is strongest in electrically neutral 
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bilayers i ' i ' Initial studies carried out at the Hartree- 
Fock level 7 identified a family of competing pseudospin 
ferromagnet states. Perturbative renormalization group 
calculation s 8,21 ' 22 have confirmed that these instabilities can 
survive beyond mean-field-theory. 
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FIG. 2: Left panel: Non-interacting biased tight-binding model 
(dotted line) and unbiased Hartree-Fock (solid line) band structures 
for bilayer graphene for k moving away from the T point towards 
the Dirac point K. The broken symmetry state bands are compared 
with non-interacting electron bands with inversion symmetry explic- 
itly broken by an externally applied electric field £ = OAV/nm. 
The interacting system bands exhibit enhanced velocities and a band 
gap due to spontaneously broken inversion symmetry. These results 
were obtained for a model with hopping parameters "fy = — 3.12eV, 
7i = — 0.377eV, y$ = j\ = 0, and dielectric screening parameter 
£ r = 4. Right panel: Onsite ^-dependent exchange potentials of the 
spontaneously broken symmetry state on the four bilayer sublattices. 
The onsite potential is larger in magnitude on the low-energy sites 
that do not have an opposite layer neighbor, and lower on average in 
the layer with the larger density. 



Fig. [2] illustrates typical mean-field theory band struc- 
tures for a 7r-band tight-binding model of unbiased bilayer 
graphene. The non-interacting bands exhibit the p 2 dispersion 
at small p which is captured by the massive chiral fermion 
model. When interaction effects are included in the mean- 
field Hamiltonian, a gap opens up and velocities increase sub- 
stantially. The gap to the remote bands associated with the 
high energy states is also increased. A gap at p = (where 
k = K so that /(k) vanishes) is always associated with a dif- 
ference in mean-field atomic ^-orbital energies between the 
two low-energy sites A and B and therefore a violation of 
the inversion symmetry which makes the two layers equiva- 
lent. We refer to this broken symmetry state as a pseudospin 
ferromagnet, motivated by its continuum model description 7 
in which the sublattice degree-of-freedom plays the role of 
a pseudospin. The increase in velocity at non-zero p illus- 
trated in Fig. [2] occurs even when inversion symmetry is not 
broken 25 . In the same figure we plot the ^-dependent on- 
site exchange potentials of the broken symmetry state which 



have been obtained from the self consistent solutions. The 
exchange potentials have opposite signs on opposite layers, 
larger magnitude at the low energy sites A and B, and values 
that increase as the valley points are approached. 

It is appropriate at this point to address some of the diffi- 
culties which arise in constructing reliably predictive theories 
for the influence of electron-electron interactions in graphene 
sheets. In continuum model theories it is customary to in- 
troduce a relative dielectric constant e r ~ (£ su h + l)/2 to ac- 
count for dielectric screening due to the substrate on which the 
graphene sheet lies. (e su b is the dielectric constant of the sub- 
strate. e su i, ~ 4 for substrates commonly used to support me- 
chanically exfoliated graphene samples. This effect is impor- 
tant for some graphene sheet properties but is often omitted in 
ab initio calculations.) Some of the results we present below 
suggest that the continuum model is reliable, even for neutral 
graphene systems (which have less screening) and even when 
E r is small, as it should be in a model intended to describe 
suspended graphene samples and in models of graphene on a 
substrate with a small dielectric constant. As we explain later, 
the two-band continuum model tends to be less accurate for 
bilayer graphene than for single-layer graphene. The impli- 
cation for the lattice model employed here, is that the on-site 
Coulomb interaction, determined by the charge form factor, 
will have a bearing on the model's predictions especially when 
S r is small. Hartree-Fock mean-field theory, used in this pa- 
per to capture the non-local exchange properties which drive 
pseudospin ferromagnetism, tends to overestimate the onset of 
broken symmetries. In ab initio calculations it is common^ to 
reduce the strength of bare exchange, say by factor of ~ 2, to 
account for Coulomb correlation screening missing in an ex- 
change only theory. This type of consideration may justify us- 
ing a value of e r larger than the one which would be suggested 
by dielectric screening considerations alone, adding another 
level of uncertainty to any quantitative predictions. 



A. Total layer density and Chern number classifications of 
competing states 

Because each bilayer flavor can polarize toward either of 
the two layers, there are a total number of 2 4 = 16 possible 
configurations of the broken symmetry 26 state. The layer po- 
larization for each spin and valley determines the sign of the 
mass term in its continuum model and has implications for the 
topological properties in the system as we will now discuss. 
The 16 states can be classified 7 by overall layer polarization 
as being ferromagnetic, ferrimagnetic, or layer antiferromag- 
netic. By this classification 7 there are two layer ferromagnetic 
states in which all flavors choose the same polarization, eight 
layer ferrimagnetic states in which three of the four flavors 
choose the same polarization, and six layer antiferromagnetic 
states with no overall polarization. 10 The layer antiferromag- 
netic states are electrostatically favored in the absence of a 
potential bias. 

One of the most interesting properties of gapped bilayer 
graphene is the existence of a finite Hall conductivity and or- 
bital magnetism due to the flavor-dependent momentum-space 
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TABLE I: Summary of Hall transport properties for valley T z = K,K' and spin G z =t,4- dependent layer polarization At = T,B. These states 
are classified as layer ferromagnetic (F), layer ferrimagnetic (Fi) or layer antiferromagnetic (AF) following reference I7I1CI1 . F, Fi, and AF 
states have respectively four, three, and two valley-spin components polarized toward the same layer (A, = T). We have listed only eight out 
of a total of sixteen configurations omitting the equivalent configurations obtained by layer reversal for every flavor. Each valley contributes a 
finite Hall conductivity with magnitude Kr and a sign that reverses with both valley interchange and layer polarization. The Hall conductivities 
are assigned to particular valleys on the basis of approximate Chern indices obtained by integrating Berry curvatures over the portion of the 
Brillouin zone near K or K' . The edge-state structure is expected to depend on each of the partial Hall conductivities. The total Hall conductivity 
is separated into contributions from separate spins and separate valleys using rj(.°' = rjj v + G^y — + ®fy - The F configuration have zero 
total hall conductivity, whereas all Fi configurations have a finite total anomalous Hall conductivity of two units. The A configurations are 
expected to be electrostatically favored in the absence of an external layer bias potential and include three types of solutions with different Hall 
properties, one with a total anomalous total Hall conductivity, one with zero Hall conductivity but finite spin Hall conductivity, and another 
with zero Hall and spin Hall conductivities because of opposite Hall conductivity contributions from K and K' valleys for both spins. 
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in the broken symmetry states. Because the vor- 
ticity v is opposite for opposite valleys, the integrated Berry 
curvature gives rise to a Hall conductivity 11 - 12 - 28 - 29 with mag- 
nitude e 2 /h for each flavor, and a sign that changes with val- 
ley as well as with layer polarization. The Berry curvature 
reflects the handedness of Bloch electrons and captures intra- 
cell circulating currents which generate a finite orbital mag- 
netic moment proportional to the angular momentum due to 
self-rotating Bloch wave packets. 29 The Berry curvature of the 
system can be evaluated usingii 
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where \u„) represents the Bloch eigenstates of the system and 
E n are the associated eigenvalues for each k. The finite or- 
bital moment generated by these wave packets has a similar 
expression and can be evaluated through 11 
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A weak external magnetic field will tend to favor a state 
in which the orbital magnetizations of all valleys are aligned, 
and the anomalous Hall conductivity is correspondingly max- 
imized. This observation suggests the possibility of valley 
optoelectronics that exploits the circular dichroism of inter- 
band transitions— in the broken symmetry states. The Kerr 
and Faraday effect measurements with linearly polarized light 



can be a useful tool to detect signatures of broken time rever- 
sal symmetry. 44 In Fig. [3] we present the ^-dependent mag- 
netization evaluated for one self-consistent gapped state in 
the presence of an interlayer bias. Similar results have been 
obtained previously using the massive Dirac-fermion contin- 
uum model. An estimate of the zero field magnetization 



per valley-spin degree of freedom M T , a , = / m„(k) 
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2^Ln .f m„ (k) k y^tt can be obtained integrating the orbital 

moment around each one of the valleys for a given spin com- 
ponent. For the broken symmetry state with e, = 4 each com- 
ponent integrates to M x _ a _ = 10~ 3 /ig per carbon atom, com- 
pared to M x _ a _ = 1.4 • 10~ 3 /Ib for the non-interacting system 
with bias <S = 0.1 V /nm andM T .<j_ = 2.2 • 10~ 3 /Xb per carbon 
atom for <S = OAV/nm. The individual flavor orbital mag- 
netizations and anomalous Hall contributions cancel in states 
that do not have broken time-reversal symmetry. 

In the non-interacting electron biased bilayer graphene 
state, opposite contributions from the two valleys lead to van- 
ishing total Hall conductivity, consistent with the absence of 
broken time-reversal symmetry. The nature of the layer po- 
larized broken symmetry we describe here is analogous to the 
biased bilayer in the sense that the charge transfer can be at- 
tributed to the (£-point dependent) exchange potential differ- 
ence between low-energy sites on opposite layers, as illus- 
trated in Fig. [2] Because of the non-locality of the exchange 
interactions, however, we are able to find self-consistent so- 
lutions in which the effective interlayer bias potential has op- 
posite signs in the two valleys, and the total Hall conductivity 
is finite. In Table I we present a list of the 16 different con- 
figurations for which we have found self-consistent solutions, 
which are characterized by the sense of layer (A, = T,B) po- 
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FIG. 3: Le/f panel: Berry curvatures associated with the two valence 
bands for a Hartree-Fock (HF) solution, and for non-interacting bi- 
layer graphene in the presence of an inteiiayer bias. We use a thicker 
line to represent the curvature of the low energy band and a thinner 
line for the band farther from the Fermi level. The small remote band 
contribution has been magnified x20. An interlayer electric field of 
$ = 0.1 V jnm added to the non-interacting electron model gives a 
band gap comparable to the one which emerges from our mean field 
calculations with £,- = 4. Right panel: Orbital magnetization contri- 
butions from the two valence bands in the vicinity of a valley point. 
We compare results for the self-consistent broken symmetry states 
with those for non-interacting electrons with an external interlayer 
bias. 



larization for each valley (t 7 = K,K') and spin (av =t>-l) as 
also discussed in reference flQIl . These results suggest the in- 
teresting possibility of altering the quantum Hall conductivity 
of a graphene bilayer sample with an external bias potential as 
we discuss later. 



B. Broken inversion symmetry states in an unbiased bilayer 

We start by examining the layer antiferromagnetic charge 
balanced configurations which are lowest in energy in the ab- 
sence of an external bias because 7 of the absence of a Hartree 
energy penalty. As summarized in Table I, there are three dis- 
tinct types of layer antiferromagnets. 10 One configuration is 
the anomalous Hall (AH) state, with four quantized units of 
Hall conductivity, in which electrons are polarized towards 
the same layer for both spins 9 , and towards opposite lay- 
ers for opposite valleys. Second is the spin Hall (SH) state 
which has opposite layer polarization on opposite valleys, and 
in each valley opposite layer polarization for opposite spin. 
The mean-field Hamiltonian for this state is similar to that of 
a non-interacting system with intrinsic spin orbit coupling 31 . 
Finally there is a solution with the same layer polarization in 
the two valleys, but opposite layer polarizations for opposite 
spins. This state has a finite valley Hall (VH) conductivity 
for each spin, zero valley Hall conductivity when summed 



over spins, and zero total Hall conductivit y '' . Optical mea- 
surements based on linearly polarized light using the Kerr and 
Faraday effect can be used to identity the anomalous Hall state 
which has broken time reversal symmetry stated 

In the continuum model formulation the three distinct solu- 
tions are degenerate. In the lattice model Hartree-Fock theory 
the energies of the AH and SH states are still exactly degen- 
erate because the exchange energy is spin diagonal and the 
energy within each spin is independent of the Hall effect sign. 
We refer to these two states collectively as the anomalous Hall 
states. The valley Hall (VH) solution does have a different en- 
ergy however, because the relative sense of layer polarization 
in the two valleys influences inter-valley exchange potentials. 
Table [II] presents our results for the condensation energy of 
the broken symmetry state and for the differences in energy 
between the VH and AH/SH solutions separated into differ- 
ent contributions. We first note that the condensation energies 
are reasonably independent of the model parameters e r and 
ro, which are not precisely known and dependent on the sam- 
ple's dielectric environment. The scale of the condensation 
energy should be compared with the value of the Coulomb in- 
teraction at the momentum scale k* (e 2 k*) over which the two- 
band model applies to bilayer graphene. We find k* by setting 
the low-energy model parabolic band energy equal to the iso- 
lated layer energy, which gives k* ~ J\/Tiv. The Coulomb en- 
ergy at this momentum scale is ~ (e 2 /e r hv) x j\ = a gr x yi ~ 
lOOmeV. This energy scale gives an estimate of the size of 
the gap which would be expected if all states with k < k* 
were fully layer polarized. We find a gap which is approx- 
imately four times smaller and density shifts between layers 
that are smaller than nk* 2 /(2n) 2 ~ (j\/2nWa) 2 by a simi- 
lar fraction. These numerical values indicate that the broken 
symmetry arises mainly from those bilayer band states that are 
reasonably well described by the low-energy two-band mod- 
els used in the original mean-field calculations, as expected. 
We note that the condensation energy is smaller by several 
orders of magnitude than the total interaction energy, which 
involves electrons across the full 7T-band. 

The second block of columns present differences in energy 
between VH states and AH/SH states. For the parameters we 
have considered the total energies for the VH states are lower 
than the AH/SH thanks to the exchange energy advantage of 
the former states. The differences in energy between these 
states become substantially smaller than the condensation en- 
ergy when the carbon radius model parameter ro is increased. 
Larger values of ro improve the accuracy of the continuum 
model and differences in energy between the different solu- 
tions rapidly decrease. We have previously argued that values 
of ro ~ 3«o m the form factor are appropriate. If so, the differ- 
ences in energy between different solutions are ~ l(T 9 eV per 
carbon atom, about 100 times smaller than the ordering con- 
densation energy. We have also explored an alternative formu- 
lation of the Hartree-Fock equations which evaluates the non- 
local exchange potentials in real space. This formulation al- 
lows the model's onsite repulsion U to be adjusted separately 
from the longer ranged tail a 33 ' 34 and thus allow a more direct 
assessment of the impact of lattice scale details of the effective 
Coulomb interaction, and leads to similar conclusions. 
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TABLE II: Energy differences between distinct mean-field solutions in the absence of a potential bias. The first group of columns represent 
energy differences between the broken symmetry pseudospin ferromagnet states in the anomalous Hall (AH) configuration and the unbroken 
symmetry state that has no gap. The latter has been obtained constraining the self-consistent HF calculation to preserve inversion symmetry. 
The third column represents the exchange energy difference AE% B of the broken symmetry self-consistent state with respect to the reference 
exchange energy of the non-interacting state. The remaining columns represent energy differences between between anomalous Hall and 
valley Hall (VH) solutions. The calculated energy differences depend sensitively on the choice of the model parameter ro used in the form 
factor. The total exchange energy differences AE'-^' = 4(AE^ K + AE% K ) can be expressed as a sum of intravalley (KK) and intervalley (KK') 
contributions, where the four-fold factor is due to the twofold degeneracy in both spin and valley. We have verified that the energy differences 
depend very weakly on the 73 and 74 parameters which are excluded in the minimal model. The two anomalous Hall states have the same 
energy in mean-field theory as explained in the text. The values compiled in this Table have been evaluated using £ r = 4 and the indicated 
carbon atom radii ro. All energies are expressed in eV per carbon atom unit. 





E (AH/SH)_ Eq 


E (VH) _ £ (AH/SH) 


ro 


AE fo! 


AE X 


A£™ 


&Etot 


AE X 


AE* K 




ao 
2a 
3ao 


-4.84- irr 8 
-5.17- ltr 8 

-6.17- 1(T 8 


-1.13- ltr 7 

-1.24- 10- 7 

-1.58- ltr 7 


-5.10- irr 4 
-6.15 -ltr 4 
-5.09 -ltr 4 


-1.03- 10 7 
-4.66- 10~ 9 
— 1.14 • 10^ 9 


-9.72- 10~ 7 
-1.92- 10~ 8 
-2.76- 10~ 9 


-1.75 -lO" 7 
-4.09- 10~ 9 
-5.71 -lO" 10 


-6.75 -10" 8 
-6.97- 10~ 10 
1.18 - 10 10 



TABLE III: Transferred charge per valley-spin flavor in the charge balanced anomalous Hall and valley Hall solutions. In our mean field 
Hartree-Fock calculations the VH solutions have slightly larger band gaps and interlayer charge transfers than their anomalous Hall counter- 
parts. In this table, charge densities are in units of 10 11 cm~ 2 and band gaps in meV units. We have used the tight-binding model parameters 
70 = — 3.l2eV, 71 = — 0.377 eV, 73 = — 0.29eV, and 74 = — 0A2eV . The 73 parameter has only a marginal influence on the broken symmetry 
state. The value chosen for 74 term tends to accumulate more electrons at the low energy sites in the bilayer but does not influence the strength 
of the broken symmetry. The results reported here were obtained using dielectric constant e r = 4 and carbon atom radius ro = 3oo- In the 
rightmost set of columns we show the results obtained when a smaller carbon atom radius ro = 2«o is used in the form factor calculation in 
Eq. ©. 
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0.46 
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0.34 


5.20 


-4.86 


-4.63 


4.29 


16 
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0.42 
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Valley Hall (VH) 
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Large values of U imply strong short-range correlations 
which are not accurately described by the continuum model. 
From the results shown in table [II] however, we see that only 
moderate short-range screening from degrees-of-freedom out- 
side the 7£-band model are required to make deviations from 
the continuum model small. Given the abundance of separate 
evidence that interaction effects in graphene systems are ac- 
curately described by a continuum model, we assume that the 
required short-range screening is in fact present. 

In Table [HI] we report on the layer distribution of charge 
when inversion symmetry breaking occurs within a particular 
valley. The results presented here are single valley-spin re- 



sults based on the AH and VH solutions. We have performed 
these calculations for band models which include and exclude 
tight-binding hopping parameters other than in plane nearest 
neighbor hopping 70 and interlayer tunneling yi ■ The trigonal 
warping 73 term connects the sites A and B and introduces a 
triangular distortion in the band structure near the Dirac point, 
hence breaking the approximately circular symmetry of the 
bands. It also makes the dispersion become linear instead of 
quadratic at the lowest energy scales and works against sym- 
metry breaking. We notice that the reduction of the band gap 
due to trigonal warping varies between two per cent and ten 
per cent depending on the short-range interaction strength, in- 
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dicating that trigonal warping plays a relatively minor role. 
The 74 term increases the accumulation of charge density at 
lattice sites A and B relative to that at the high energy sites 
B and B. The gap size and total transferred charge (An/) re- 
main virtually unchanged, although the distribution of charge 
between sublattices on the same layer (An l s ) is altered in some 
states. 
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FIG. 4: Temperature dependence of mean-field-theory charge trans- 
fer per valley-spin and the associated band gap. Top panel: Total 
charge per unit area (10 11 cm~ 2 ) transferred from one layer to an- 
other An; per each valley-spin polarized component as a function of 
temperature. We represent the dependence for different values of 
the relative dielectric constant e r . Middle panel: The bottom panel 
presents the mean-field-theory band gap as a function of temperature. 
Bottom panel: Ratio between the band gap and transferred charge 
density. 



The temperature dependence of the band gap and the trans- 
ferred charge per valley-spin component in mean-field-theory 
is plotted in Fig. [4] As the temperature is increased both 
the charge density and band gaps are reduced, but their ra- 
tio is approximately fixed. The decay trend is similar for the 
different dielectric constants we have considered. The criti- 
cal temperatures obtained in this mean-field calculation pro- 
vide an estimate of the maximum temperature at which local 
order survives in the absence of disorder. The band gap de- 
creases more quickly than the charge-density transfer order 
parameter as the interaction is made weaker; the charge den- 
sity per valley-spin varies approximately like e,r ! whereas the 
gap varies as e ( 7 155 . 



Quantum phase transitions in the low magnetic field and 
low bias regime 



To first order in magnetic field, the magnetic contribution 
to energy is simply proportional to spontaneous magnetization 
discussed previously, which is entirely orbital in character and 
depends on the flavor-dependent layer polarizations. States 
that are related by layer polarization reversal are no longer 
equal in energy in the presence of a magnetic field. If we 
take AE tot ~ 10~ 9 eV per carbon atom for the mean field en- 
ergy differences between AH and VH states from TableHIland 
use the relation AE tot ~ M T _ a7 ■ B c , we obtain from the orbital 
magnetization values we have calculated that a magnetic field 
B c ~ 0.004r is sufficient to favor the AH state with orbital 
magnetization parallel to the magnetic field over VH states. 
Considering the small energy differences between the differ- 
ent competing states as shown in table [II] we can expect that 
the coupling between magnetic field and the orbital moment in 
the system can play a decisive role in selecting the minimum 
energy ground-state. Because the electron densities at which 
gaps occur are magnetic field dependent in states with finite 
Hall conductivity, small fluctuations in the density like those 
associated with electron-hole puddle domains 42 will also have 
an important influence on the nature of the ground-state at low 
magnetic fields. 

We now consider the case of zero magnetic field and finite 
electric field. In presence of an external bias the inversion 
symmetry of a graphene bilayer is explicitly broken, favoring 
charge accumulation in one of the layers. When the exter- 
nal bias becomes large enough the charge balanced broken 
symmetry antiferromagnetic (AF) configuration gives way to 
ferrimagnetic (Fi) and eventually ferromagnetic (F) configu- 
rations. These solutions can become energetically favored for 
certain ranges of the external bias potential thanks to their in- 
trinsic broken inversion symmetry configuration with sponta- 
neous charge transfer. In Fig. [3] we present the total ener- 
gies of the different types of solutions in the low bias regime, 
obtained by starting the self-consistent calculations from dif- 
ferent seeds. For low potential bias the charge balanced AF 
structure remains lowest in energy. Eventually the external 
field becomes large enough to flip the layer polarization of one 
valley-spin component towards the layer favored by the elec- 
tric field, giving rise to the Fi type solutions. The electrostatic 
charge imbalance of the Fi states is approximately half of that 
associated with the F state in which all four components are 
polarized towards the same layer, a fact that can be inferred 
from the slopes of the total energy evolution as a function of 
external electric field. In Fig. [6] we present the charge densi- 
ties associated with one valley-spin component with different 
layer polarizations. For the AF and Fi solutions we have com- 
ponents that are polarized toward both layers that we represent 
as An/ = («; — no), the density difference with respect to the 
uniform background density nrj. In the Fi configuration the 
three charge density components polarized towards the same 
layer do not have exactly the same value when the electron 
spins are different, but they are similar in magnitude. In the 
F configuration all four density components are polarized to- 
wards the same layer and have the same magnitude. 
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FIG. 5: Bias dependence of the total energies of the system for 
different valley-spin polarization configurations for a system with 
^ r = 4. The system undergoes transitions from the AF to Fi and then 
to F states as a function of external field, each one displaying clearly 
different quantum Hall conductivity properties. The origin of energy 
has been arbitrarily shifted for presentation convenience. 



The behavior of the band gaps in the low bias regime illus- 
trates the character of each solution. The AF band gap con- 
tinually decreases when the electric field is increased. In this 
case the electric field always works against the charge density 
distribution of two components directly contributing in the re- 
duction of the band gap. For the F configuration we have a 
completely opposite trend; the gap always increases since the 
electric field favors charge accumulation in the layer in which 
the density is already higher. In the Fi configuration we find 
an intermediate situation. We find an initial increase of the 
gap in presence of an electric field thanks to the reduction 
of the Hartree energy penalty associated to the spontaneous 
charge imbalance. Beyond a certain point the band gap starts 
to decrease, reflecting the fact that the electric field is working 
against the layer polarization of one component. Eventually 
when the effect of the external field becomes dominant the 
system undergoes a phase transition to the F configuration. 
Every time there is a crossover of the minimum energy lev- 
els we will see an abrupt change in the band gap and charge 
density. The Hall conductivity of the system will also see 
discontinuous changes following the classification of table U 
The observations above suggest that the band gap will show 
a non-monotonic dependence on potential bias with an initial 
reduction at low fields for AF states and an increase for large 
enough bias for F states. Between those two limits the system 
can form a Fi state that also has a band gap that decreases with 
increasing electric field, before finally making a transition to 
the F state. 



D. Exchange screening effects in presence of a strong bias 

In the presence of an external bias the inversion symmetry 
of graphene bilayer is explicitly broken, favoring the charge 
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FIG. 6: Upper panel: Each branch represents the bias dependence 
of the valley resolved layer polarization density associated with one 
valley-spin component. These results were obtained for a system 
with £,. = 4. The multiplication factors represent the number of times 
each branch needs to be repeated to give the total charge density to 
complete the contributions from the four components. Lower panel: 
Band gap as a function of bias for F, Fi and AF states for a sys- 
tem with e r = 4. We observe a competition in the band-gap opening 
between exchange and Hartree contributions which can lead to a re- 
duction of the band gap as the bias is increased. 



accumulation in one of the layers and opening a band gap, a 
fact that has been verified in several experiments 35 " 37 and also 
predicted theoretically within tight-binding 38 , tight-binding 
plus Hartree screenin g 39 ' 40 or ab-initio calculations^. Here 
we explore the role that electron exchange can play when the 
system is subject to a strong external bias. In Fig. [7]we present 
the charge accumulation in one of the layers as a function of 
bias. In the strong bias limit we are in the F configuration 
in which all four components are polarized towards the same 
layer. Results obtained in the Hartree approximation follows 
a similar trend in the strong bias limit resulting in compara- 
ble amounts of total transferred charge in the range of bias 
we considered, with the Hartree only screening allowing more 
sloshed charge. A larger value of dielectric screening e, that 
weakens the interaction strength also weakens the electrostatic 
screening and therefore more charge imbalance per layer is 
expected. 

As we show in Fig. |7]the presence of the exchange term in 
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the Hamiltonian introduces a clear enhancement of the band 
gap that persists up to high bias potentials. This gap enhance- 
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FIG. 7: External bias dependence of layer resolved charge density, 
band gaps and sublattice charge asymmetry. We have considered the 
F configuration where all four valley-spin components are polarized 
towards the same layer. Upper panel: Charge accumulation depen- 
dence as a function of bias in presence of Hartree and Fock terms 
and charge accumulation dependence as a function of bias in pres- 
ence of Hartree screening only. The evolution of An/ as a function 
of bias has a slightly larger slope in the Hartree only approximation 
than the Hartree-Fock case. Middle panel: Band gaps obtained for 
Hartree-Fock and Hartree approximations using different values of 
£,• as a function of external electric field S. We notice a substantial 
enhancement of the band gap when the exchange term is included. 
Lower panel: Measure of deviation from neutrality in the charge 
differences for each one of the sublattices in a given layer normal- 
ized by the total sloshed charge. The asymmetric distribution of the 
charge between low energy sublattice A and high energy sublattice 
B of a given graphene layer is enhanced in presence of electron ex- 
change. This asymmetry becomes smaller as the external bias be- 
comes stronger. 

ment can be related to the exchange contribution that tends 
to introduce a strong asymmetry in charge distribution be- 
tween A and B sublattices of a given layer that persists up 
to rather high values of external electric field as shown in Fig. 
Q We may define this intralayer charge imbalance ratio as 
(|Ara^| + |Anf |) /|A«/|. The enhancement of charge imbal- 
ance between the sublattices is one of the effects of intralayer 
exchange that contributes to changing the magnitude of the 
band gap. This quantity is more sensitive to the strength of 
the Coulomb interaction than to changes in the temperature. 



V. SUMMARY AND CONCLUSIONS 

We have presented a numerical description of the broken 
charge symmetry Hartree-Fock ground state in a graphene bi- 
layer system for different dielectric screening of the media 
and temperature. The predicted gap sizes and temperature 
ranges in which this effect is expected to survive should be 
experimentally accessible. The size of the interlayer charge 
transfers suggests that the broken symmetry states will be sup- 
pressed by smooth disorder strong enough to produce charge 
puddles with density variations larger than around 10~ 5 per 
carbon atom. The total amount of spontaneous charge trans- 
fer from one layer to another in the broken ground state is 
~ 10 1 l cm~ 2 , smaller but comparable in magnitude to the den- 
sities that can be induced by gating the device or depositing 
impurities 1 and can therefore be relevant for interpreting gat- 
ing experiments in bilayer graphene. The mechanism for bro- 
ken symmetry described here is most effective for a neutral bi- 
layer because the energy gain due to charge transfer is greatest 
for electrons nearest to the Dirac point. The broken symme- 
tries are expected to be present only when the carrier density is 
smaller than about 10 11 cm -2 , hence a reduction of this effect 
is expected as the Fermi surface moves away from the charge 
neutrality point when the system is doped. 

The Coulomb correlation neglected in the Hartree-Fock for- 
malism usually plays a more important in extended systems 41 
than in localized atomic systems where they are typically one 
order of magnitude smaller than exchange effects. Treatments 
of Coulomb correlation based on a renormalization group 
framework 8 or RPA 23 suggest that the broken symmetry ad- 
dressed in this paper will survive when Coulomb correlation 
is considered. In our case we have used larger values of S r 
than those based on estimates from dielectric functions alone 
in order to effectively reduce the strength of exchange to ap- 
proximately account for screening effects. Application of a 
self-consistent dynamic screening approximation to this prob- 
lem would be interesting but would require the high density 
of sampling ^-points near the Dirac point to be implemented 
in a computationally efficient way. 

An interesting feature of these broken symmetry states is 
the spontaneous quantum Hall effects which appear when op- 
posite valleys have opposite layer polarizations^ When op- 
posite valleys have the same layer polarization, the system has 
only a valley Hall effect, which is not manifested in standard 
electrical measurements. There is no energy difference be- 
tween these states in a continuum model which does not ac- 
count for lattice-scale physics. We find that the energy differ- 
ence between these states decreases rapidly depending on the 
strength of the on-site effective interaction in our 7T-band only 
model, with valley Hall states being favored over anomalous 
Hall states. For a very weak on site interaction, anomalous 
Hall states could be favored over valley Hall states. States in 
which electrons with one spin orientation form one Hall state, 
while electrons with the other spin form the other Hall state 
also occur. Despite the uncertainties of the calculation due to 
this high sensitivity to lattice details they do give us an es- 
timate for the magnitude of the energy differences between 
different configurations. Unlike Nandkishore and Levitov, 9 
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we find that valley Hall states have lower energy than anoma- 
lous Hall states, but that the energy competition can be altered 
by even quite weak external magnetic fields. This scenario 
appears to be consistent with experiments from the Yacoby 
group in which a v = 4 quantized Hall effect persists down 
to very weak magnetic fields. The small energy differences 
between competing states could mean that domains with all 
characters are present, separated by domain walls, because 
of entropic and disorder considerations. In the presence of 
a magnetic field, coupling of spontaneous orbital magnetism 
to the external field favors anomalous Hall states and should 
coarsen any domain structure. Similarly an externally applied 
potential bias favors valley Hall states and should also coarsen 
domain structures. In the absence of disorder we do find that 
the charge gap is first decreased by potential bias, and finally 
increased once both spins have valley Hall effects with the 
same layer polarization. This finding is qualitatively consis- 
tent with experiment^ 

We have also presented mean field theory estimates of the 



critical temperatures associated with spontaneous layer polar- 
ization states in graphene bilayers. According to these esti- 
mates, order will not survive to room temperatures and there- 
fore is unlikely to be useful for applications. On the other 
hand, the exchange interactions which produce broken sym- 
metry states at low temperatures, will still play a role in en- 
hancing gaps produced by external potential biases at room 
temperature. 
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